home *** CD-ROM | disk | FTP | other *** search
/ Mac Format 1995 June / MacFormat 25.iso / Shareware City / Developers / fortran-to-c-translator-11 / Mac F2C 1.1 / Mac F2C Documentation / If Floats Don't Display Right / dtoa.c next >
Text File  |  1994-12-18  |  48KB  |  2,475 lines
  1. /****************************************************************
  2.  *
  3.  * The author of this software is David M. Gay.
  4.  *
  5.  * Copyright (c) 1991 by AT&T.
  6.  *
  7.  * Permission to use, copy, modify, and distribute this software for any
  8.  * purpose without fee is hereby granted, provided that this entire notice
  9.  * is included in all copies of any software which is or includes a copy
  10.  * or modification of this software and in all copies of the supporting
  11.  * documentation for such software.
  12.  *
  13.  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
  14.  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
  15.  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
  16.  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
  17.  *
  18.  ***************************************************************/
  19.  
  20. /* Please send bug reports to
  21.     David M. Gay
  22.     AT&T Bell Laboratories, Room 2C-463
  23.     600 Mountain Avenue
  24.     Murray Hill, NJ 07974-2070
  25.     U.S.A.
  26.     dmg@research.att.com or research!dmg
  27.  */
  28.  
  29. #define IEEE_MC68k
  30.  
  31. /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
  32.  *
  33.  * This strtod returns a nearest machine number to the input decimal
  34.  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
  35.  * broken by the IEEE round-even rule.  Otherwise ties are broken by
  36.  * biased rounding (add half and chop).
  37.  *
  38.  * Inspired loosely by William D. Clinger's paper "How to Read Floating
  39.  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
  40.  *
  41.  * Modifications:
  42.  *
  43.  *    1. We only require IEEE, IBM, or VAX double-precision
  44.  *        arithmetic (not IEEE double-extended).
  45.  *    2. We get by with floating-point arithmetic in a case that
  46.  *        Clinger missed -- when we're computing d * 10^n
  47.  *        for a small integer d and the integer n is not too
  48.  *        much larger than 22 (the maximum integer k for which
  49.  *        we can represent 10^k exactly), we may be able to
  50.  *        compute (d*10^k) * 10^(e-k) with just one roundoff.
  51.  *    3. Rather than a bit-at-a-time adjustment of the binary
  52.  *        result in the hard case, we use floating-point
  53.  *        arithmetic to determine the adjustment to within
  54.  *        one bit; only in really hard cases do we need to
  55.  *        compute a second residual.
  56.  *    4. Because of 3., we don't need a large table of powers of 10
  57.  *        for ten-to-e (just some small tables, e.g. of 10^k
  58.  *        for 0 <= k <= 22).
  59.  */
  60.  
  61. /*
  62.  * #define IEEE_8087 for IEEE-arithmetic machines where the least
  63.  *    significant byte has the lowest address.
  64.  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
  65.  *    significant byte has the lowest address.
  66.  * #define Long int on machines with 32-bit ints and 64-bit longs.
  67.  * #define Sudden_Underflow for IEEE-format machines without gradual
  68.  *    underflow (i.e., that flush to zero on underflow).
  69.  * #define IBM for IBM mainframe-style floating-point arithmetic.
  70.  * #define VAX for VAX-style floating-point arithmetic.
  71.  * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
  72.  * #define No_leftright to omit left-right logic in fast floating-point
  73.  *    computation of dtoa.
  74.  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
  75.  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
  76.  *    that use extended-precision instructions to compute rounded
  77.  *    products and quotients) with IBM.
  78.  * #define ROUND_BIASED for IEEE-format with biased rounding.
  79.  * #define Inaccurate_Divide for IEEE-format with correctly rounded
  80.  *    products but inaccurate quotients, e.g., for Intel i860.
  81.  * #define Just_16 to store 16 bits per 32-bit Long when doing high-precision
  82.  *    integer arithmetic.  Whether this speeds things up or slows things
  83.  *    down depends on the machine and the number being converted.
  84.  * #define KR_headers for old-style C function headers.
  85.  * #define Bad_float_h if your system lacks a float.h or if it does not
  86.  *    define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
  87.  *    FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
  88.  * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
  89.  *    if memory is available and otherwise does something you deem
  90.  *    appropriate.  If MALLOC is undefined, malloc will be invoked
  91.  *    directly -- and assumed always to succeed.
  92.  */
  93.  
  94. #ifndef Long
  95. #define Long long
  96. #endif
  97.  
  98. #ifdef DEBUG
  99. #include "stdio.h"
  100. #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
  101. #endif
  102.  
  103. #ifdef __cplusplus
  104. #include "malloc.h"
  105. #include "memory.h"
  106. #else
  107. #ifndef KR_headers
  108. #include "stdlib.h"
  109. #include "string.h"
  110. #else
  111. #include "malloc.h"
  112. #include "memory.h"
  113. #endif
  114. #endif
  115.  
  116. #ifdef MALLOC
  117. #ifdef KR_headers
  118. extern char *MALLOC();
  119. #else
  120. extern void *MALLOC(size_t);
  121. #endif
  122. #else
  123. #define MALLOC malloc
  124. #endif
  125.  
  126. #include "errno.h"
  127. #ifdef Bad_float_h
  128. #undef __STDC__
  129. #ifdef IEEE_MC68k
  130. #define IEEE_ARITHMETIC
  131. #endif
  132. #ifdef IEEE_8087
  133. #define IEEE_ARITHMETIC
  134. #endif
  135.  
  136. #ifdef IEEE_ARITHMETIC
  137. #define DBL_DIG 15
  138. #define DBL_MAX_10_EXP 308
  139. #define DBL_MAX_EXP 1024
  140. #define FLT_RADIX 2
  141. #define FLT_ROUNDS 1
  142. #define DBL_MAX 1.7976931348623157e+308
  143. #endif
  144.  
  145. #ifdef IBM
  146. #define DBL_DIG 16
  147. #define DBL_MAX_10_EXP 75
  148. #define DBL_MAX_EXP 63
  149. #define FLT_RADIX 16
  150. #define FLT_ROUNDS 0
  151. #define DBL_MAX 7.2370055773322621e+75
  152. #endif
  153.  
  154. #ifdef VAX
  155. #define DBL_DIG 16
  156. #define DBL_MAX_10_EXP 38
  157. #define DBL_MAX_EXP 127
  158. #define FLT_RADIX 2
  159. #define FLT_ROUNDS 1
  160. #define DBL_MAX 1.7014118346046923e+38
  161. #endif
  162.  
  163. #ifndef LONG_MAX
  164. #define LONG_MAX 2147483647
  165. #endif
  166. #else
  167. #include "float.h"
  168. #endif
  169. #ifndef __MATH_H__
  170. #include "math.h"
  171. #endif
  172.  
  173. #ifdef __cplusplus
  174. extern "C" {
  175. #endif
  176.  
  177. #ifndef CONST
  178. #ifdef KR_headers
  179. #define CONST /* blank */
  180. #else
  181. #define CONST const
  182. #endif
  183. #endif
  184.  
  185. #ifdef Unsigned_Shifts
  186. #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
  187. #else
  188. #define Sign_Extend(a,b) /*no-op*/
  189. #endif
  190.  
  191. #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
  192. Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
  193. #endif
  194.  
  195. #ifdef IEEE_8087
  196. #define word0(x) ((unsigned Long *)&x)[1]
  197. #define word1(x) ((unsigned Long *)&x)[0]
  198. #else
  199. #define word0(x) ((unsigned Long *)&x)[0]
  200. #define word1(x) ((unsigned Long *)&x)[1]
  201. #endif
  202.  
  203. /* The following definition of Storeinc is appropriate for MIPS processors.
  204.  * An alternative that might be better on some machines is
  205.  * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
  206.  */
  207. #if defined(IEEE_8087) + defined(VAX)
  208. #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
  209. ((unsigned short *)a)[0] = (unsigned short)c, a++)
  210. #else
  211. #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
  212. ((unsigned short *)a)[1] = (unsigned short)c, a++)
  213. #endif
  214.  
  215. /* #define P DBL_MANT_DIG */
  216. /* Ten_pmax = floor(P*log(2)/log(5)) */
  217. /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
  218. /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
  219. /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
  220.  
  221. #if defined(IEEE_8087) + defined(IEEE_MC68k)
  222. #define Exp_shift  20
  223. #define Exp_shift1 20
  224. #define Exp_msk1    0x100000
  225. #define Exp_msk11   0x100000
  226. #define Exp_mask  0x7ff00000
  227. #define P 53
  228. #define Bias 1023
  229. #define IEEE_Arith
  230. #define Emin (-1022)
  231. #define Exp_1  0x3ff00000
  232. #define Exp_11 0x3ff00000
  233. #define Ebits 11
  234. #define Frac_mask  0xfffff
  235. #define Frac_mask1 0xfffff
  236. #define Ten_pmax 22
  237. #define Bletch 0x10
  238. #define Bndry_mask  0xfffff
  239. #define Bndry_mask1 0xfffff
  240. #define LSB 1
  241. #define Sign_bit 0x80000000
  242. #define Log2P 1
  243. #define Tiny0 0
  244. #define Tiny1 1
  245. #define Quick_max 14
  246. #define Int_max 14
  247. #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
  248. #else
  249. #undef  Sudden_Underflow
  250. #define Sudden_Underflow
  251. #ifdef IBM
  252. #define Exp_shift  24
  253. #define Exp_shift1 24
  254. #define Exp_msk1   0x1000000
  255. #define Exp_msk11  0x1000000
  256. #define Exp_mask  0x7f000000
  257. #define P 14
  258. #define Bias 65
  259. #define Exp_1  0x41000000
  260. #define Exp_11 0x41000000
  261. #define Ebits 8    /* exponent has 7 bits, but 8 is the right value in b2d */
  262. #define Frac_mask  0xffffff
  263. #define Frac_mask1 0xffffff
  264. #define Bletch 4
  265. #define Ten_pmax 22
  266. #define Bndry_mask  0xefffff
  267. #define Bndry_mask1 0xffffff
  268. #define LSB 1
  269. #define Sign_bit 0x80000000
  270. #define Log2P 4
  271. #define Tiny0 0x100000
  272. #define Tiny1 0
  273. #define Quick_max 14
  274. #define Int_max 15
  275. #else /* VAX */
  276. #define Exp_shift  23
  277. #define Exp_shift1 7
  278. #define Exp_msk1    0x80
  279. #define Exp_msk11   0x800000
  280. #define Exp_mask  0x7f80
  281. #define P 56
  282. #define Bias 129
  283. #define Exp_1  0x40800000
  284. #define Exp_11 0x4080
  285. #define Ebits 8
  286. #define Frac_mask  0x7fffff
  287. #define Frac_mask1 0xffff007f
  288. #define Ten_pmax 24
  289. #define Bletch 2
  290. #define Bndry_mask  0xffff007f
  291. #define Bndry_mask1 0xffff007f
  292. #define LSB 0x10000
  293. #define Sign_bit 0x8000
  294. #define Log2P 1
  295. #define Tiny0 0x80
  296. #define Tiny1 0
  297. #define Quick_max 15
  298. #define Int_max 15
  299. #endif
  300. #endif
  301.  
  302. #ifndef IEEE_Arith
  303. #define ROUND_BIASED
  304. #endif
  305.  
  306. #ifdef RND_PRODQUOT
  307. #define rounded_product(a,b) a = rnd_prod(a, b)
  308. #define rounded_quotient(a,b) a = rnd_quot(a, b)
  309. #ifdef KR_headers
  310. extern double rnd_prod(), rnd_quot();
  311. #else
  312. extern double rnd_prod(double, double), rnd_quot(double, double);
  313. #endif
  314. #else
  315. #define rounded_product(a,b) a *= b
  316. #define rounded_quotient(a,b) a /= b
  317. #endif
  318.  
  319. #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
  320. #define Big1 0xffffffff
  321.  
  322. #ifndef Just_16
  323. /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
  324.  * This makes some inner loops simpler and sometimes saves work
  325.  * during multiplications, but it often seems to make things slightly
  326.  * slower.  Hence the default is now to store 32 bits per Long.
  327.  */
  328. #ifndef Pack_32
  329. #define Pack_32
  330. #endif
  331. #endif
  332.  
  333. #define Kmax 15
  334.  
  335. #ifdef __cplusplus
  336. extern "C" double strtod(const char *s00, char **se);
  337. extern "C" char *dtoa(double d, int mode, int ndigits,
  338.             int *decpt, int *sign, char **rve);
  339. #endif
  340.  
  341.  struct
  342. Bigint {
  343.     struct Bigint *next;
  344.     int k, maxwds, sign, wds;
  345.     unsigned Long x[1];
  346.     };
  347.  
  348.  typedef struct Bigint Bigint;
  349.  
  350.  static Bigint *freelist[Kmax+1];
  351.  
  352.  static Bigint *
  353. Balloc
  354. #ifdef KR_headers
  355.     (k) int k;
  356. #else
  357.     (int k)
  358. #endif
  359. {
  360.     int x;
  361.     Bigint *rv;
  362.  
  363.     if (rv = freelist[k]) {
  364.         freelist[k] = rv->next;
  365.         }
  366.     else {
  367.         x = 1 << k;
  368.         rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(Long));
  369.         rv->k = k;
  370.         rv->maxwds = x;
  371.         }
  372.     rv->sign = rv->wds = 0;
  373.     return rv;
  374.     }
  375.  
  376.  static void
  377. Bfree
  378. #ifdef KR_headers
  379.     (v) Bigint *v;
  380. #else
  381.     (Bigint *v)
  382. #endif
  383. {
  384.     if (v) {
  385.         v->next = freelist[v->k];
  386.         freelist[v->k] = v;
  387.         }
  388.     }
  389.  
  390. #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
  391. y->wds*sizeof(Long) + 2*sizeof(int))
  392.  
  393.  static Bigint *
  394. multadd
  395. #ifdef KR_headers
  396.     (b, m, a) Bigint *b; int m, a;
  397. #else
  398.     (Bigint *b, int m, int a)    /* multiply by m and add a */
  399. #endif
  400. {
  401.     int i, wds;
  402.     unsigned Long *x, y;
  403. #ifdef Pack_32
  404.     unsigned Long xi, z;
  405. #endif
  406.     Bigint *b1;
  407.  
  408.     wds = b->wds;
  409.     x = b->x;
  410.     i = 0;
  411.     do {
  412. #ifdef Pack_32
  413.         xi = *x;
  414.         y = (xi & 0xffff) * m + a;
  415.         z = (xi >> 16) * m + (y >> 16);
  416.         a = (int)(z >> 16);
  417.         *x++ = (z << 16) + (y & 0xffff);
  418. #else
  419.         y = *x * m + a;
  420.         a = (int)(y >> 16);
  421.         *x++ = y & 0xffff;
  422. #endif
  423.         }
  424.         while(++i < wds);
  425.     if (a) {
  426.         if (wds >= b->maxwds) {
  427.             b1 = Balloc(b->k+1);
  428.             Bcopy(b1, b);
  429.             Bfree(b);
  430.             b = b1;
  431.             }
  432.         b->x[wds++] = a;
  433.         b->wds = wds;
  434.         }
  435.     return b;
  436.     }
  437.  
  438.  static Bigint *
  439. s2b
  440. #ifdef KR_headers
  441.     (s, nd0, nd, y9) CONST char *s; int nd0, nd; unsigned Long y9;
  442. #else
  443.     (CONST char *s, int nd0, int nd, unsigned Long y9)
  444. #endif
  445. {
  446.     Bigint *b;
  447.     int i, k;
  448.     Long x, y;
  449.  
  450.     x = (nd + 8) / 9;
  451.     for(k = 0, y = 1; x > y; y <<= 1, k++) ;
  452. #ifdef Pack_32
  453.     b = Balloc(k);
  454.     b->x[0] = y9;
  455.     b->wds = 1;
  456. #else
  457.     b = Balloc(k+1);
  458.     b->x[0] = y9 & 0xffff;
  459.     b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
  460. #endif
  461.  
  462.     i = 9;
  463.     if (9 < nd0) {
  464.         s += 9;
  465.         do b = multadd(b, 10, *s++ - '0');
  466.             while(++i < nd0);
  467.         s++;
  468.         }
  469.     else
  470.         s += 10;
  471.     for(; i < nd; i++)
  472.         b = multadd(b, 10, *s++ - '0');
  473.     return b;
  474.     }
  475.  
  476.  static int
  477. hi0bits
  478. #ifdef KR_headers
  479.     (x) register unsigned Long x;
  480. #else
  481.     (register unsigned Long x)
  482. #endif
  483. {
  484.     register int k = 0;
  485.  
  486.     if (!(x & 0xffff0000)) {
  487.         k = 16;
  488.         x <<= 16;
  489.         }
  490.     if (!(x & 0xff000000)) {
  491.         k += 8;
  492.         x <<= 8;
  493.         }
  494.     if (!(x & 0xf0000000)) {
  495.         k += 4;
  496.         x <<= 4;
  497.         }
  498.     if (!(x & 0xc0000000)) {
  499.         k += 2;
  500.         x <<= 2;
  501.         }
  502.     if (!(x & 0x80000000)) {
  503.         k++;
  504.         if (!(x & 0x40000000))
  505.             return 32;
  506.         }
  507.     return k;
  508.     }
  509.  
  510.  static int
  511. lo0bits
  512. #ifdef KR_headers
  513.     (y) unsigned Long *y;
  514. #else
  515.     (unsigned Long *y)
  516. #endif
  517. {
  518.     register int k;
  519.     register unsigned Long x = *y;
  520.  
  521.     if (x & 7) {
  522.         if (x & 1)
  523.             return 0;
  524.         if (x & 2) {
  525.             *y = x >> 1;
  526.             return 1;
  527.             }
  528.         *y = x >> 2;
  529.         return 2;
  530.         }
  531.     k = 0;
  532.     if (!(x & 0xffff)) {
  533.         k = 16;
  534.         x >>= 16;
  535.         }
  536.     if (!(x & 0xff)) {
  537.         k += 8;
  538.         x >>= 8;
  539.         }
  540.     if (!(x & 0xf)) {
  541.         k += 4;
  542.         x >>= 4;
  543.         }
  544.     if (!(x & 0x3)) {
  545.         k += 2;
  546.         x >>= 2;
  547.         }
  548.     if (!(x & 1)) {
  549.         k++;
  550.         x >>= 1;
  551.         if (!x & 1)
  552.             return 32;
  553.         }
  554.     *y = x;
  555.     return k;
  556.     }
  557.  
  558.  static Bigint *
  559. i2b
  560. #ifdef KR_headers
  561.     (i) int i;
  562. #else
  563.     (int i)
  564. #endif
  565. {
  566.     Bigint *b;
  567.  
  568.     b = Balloc(1);
  569.     b->x[0] = i;
  570.     b->wds = 1;
  571.     return b;
  572.     }
  573.  
  574.  static Bigint *
  575. mult
  576. #ifdef KR_headers
  577.     (a, b) Bigint *a, *b;
  578. #else
  579.     (Bigint *a, Bigint *b)
  580. #endif
  581. {
  582.     Bigint *c;
  583.     int k, wa, wb, wc;
  584.     unsigned Long carry, y, z;
  585.     unsigned Long *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
  586. #ifdef Pack_32
  587.     unsigned Long z2;
  588. #endif
  589.  
  590.     if (a->wds < b->wds) {
  591.         c = a;
  592.         a = b;
  593.         b = c;
  594.         }
  595.     k = a->k;
  596.     wa = a->wds;
  597.     wb = b->wds;
  598.     wc = wa + wb;
  599.     if (wc > a->maxwds)
  600.         k++;
  601.     c = Balloc(k);
  602.     for(x = c->x, xa = x + wc; x < xa; x++)
  603.         *x = 0;
  604.     xa = a->x;
  605.     xae = xa + wa;
  606.     xb = b->x;
  607.     xbe = xb + wb;
  608.     xc0 = c->x;
  609. #ifdef Pack_32
  610.     for(; xb < xbe; xb++, xc0++) {
  611.         if (y = *xb & 0xffff) {
  612.             x = xa;
  613.             xc = xc0;
  614.             carry = 0;
  615.             do {
  616.                 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
  617.                 carry = z >> 16;
  618.                 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
  619.                 carry = z2 >> 16;
  620.                 Storeinc(xc, z2, z);
  621.                 }
  622.                 while(x < xae);
  623.             *xc = carry;
  624.             }
  625.         if (y = *xb >> 16) {
  626.             x = xa;
  627.             xc = xc0;
  628.             carry = 0;
  629.             z2 = *xc;
  630.             do {
  631.                 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
  632.                 carry = z >> 16;
  633.                 Storeinc(xc, z, z2);
  634.                 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
  635.                 carry = z2 >> 16;
  636.                 }
  637.                 while(x < xae);
  638.             *xc = z2;
  639.             }
  640.         }
  641. #else
  642.     for(; xb < xbe; xc0++) {
  643.         if (y = *xb++) {
  644.             x = xa;
  645.             xc = xc0;
  646.             carry = 0;
  647.             do {
  648.                 z = *x++ * y + *xc + carry;
  649.                 carry = z >> 16;
  650.                 *xc++ = z & 0xffff;
  651.                 }
  652.                 while(x < xae);
  653.             *xc = carry;
  654.             }
  655.         }
  656. #endif
  657.     for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
  658.     c->wds = wc;
  659.     return c;
  660.     }
  661.  
  662.  static Bigint *p5s;
  663.  
  664.  static Bigint *
  665. pow5mult
  666. #ifdef KR_headers
  667.     (b, k) Bigint *b; int k;
  668. #else
  669.     (Bigint *b, int k)
  670. #endif
  671. {
  672.     Bigint *b1, *p5, *p51;
  673.     int i;
  674.     static int p05[3] = { 5, 25, 125 };
  675.  
  676.     if (i = k & 3)
  677.         b = multadd(b, p05[i-1], 0);
  678.  
  679.     if (!(k >>= 2))
  680.         return b;
  681.     if (!(p5 = p5s)) {
  682.         /* first time */
  683.         p5 = p5s = i2b(625);
  684.         p5->next = 0;
  685.         }
  686.     for(;;) {
  687.         if (k & 1) {
  688.             b1 = mult(b, p5);
  689.             Bfree(b);
  690.             b = b1;
  691.             }
  692.         if (!(k >>= 1))
  693.             break;
  694.         if (!(p51 = p5->next)) {
  695.             p51 = p5->next = mult(p5,p5);
  696.             p51->next = 0;
  697.             }
  698.         p5 = p51;
  699.         }
  700.     return b;
  701.     }
  702.  
  703.  static Bigint *
  704. lshift
  705. #ifdef KR_headers
  706.     (b, k) Bigint *b; int k;
  707. #else
  708.     (Bigint *b, int k)
  709. #endif
  710. {
  711.     int i, k1, n, n1;
  712.     Bigint *b1;
  713.     unsigned Long *x, *x1, *xe, z;
  714.  
  715. #ifdef Pack_32
  716.     n = k >> 5;
  717. #else
  718.     n = k >> 4;
  719. #endif
  720.     k1 = b->k;
  721.     n1 = n + b->wds + 1;
  722.     for(i = b->maxwds; n1 > i; i <<= 1)
  723.         k1++;
  724.     b1 = Balloc(k1);
  725.     x1 = b1->x;
  726.     for(i = 0; i < n; i++)
  727.         *x1++ = 0;
  728.     x = b->x;
  729.     xe = x + b->wds;
  730. #ifdef Pack_32
  731.     if (k &= 0x1f) {
  732.         k1 = 32 - k;
  733.         z = 0;
  734.         do {
  735.             *x1++ = *x << k | z;
  736.             z = *x++ >> k1;
  737.             }
  738.             while(x < xe);
  739.         if (*x1 = z)
  740.             ++n1;
  741.         }
  742. #else
  743.     if (k &= 0xf) {
  744.         k1 = 16 - k;
  745.         z = 0;
  746.         do {
  747.             *x1++ = *x << k  & 0xffff | z;
  748.             z = *x++ >> k1;
  749.             }
  750.             while(x < xe);
  751.         if (*x1 = z)
  752.             ++n1;
  753.         }
  754. #endif
  755.     else do
  756.         *x1++ = *x++;
  757.         while(x < xe);
  758.     b1->wds = n1 - 1;
  759.     Bfree(b);
  760.     return b1;
  761.     }
  762.  
  763.  static int
  764. cmp
  765. #ifdef KR_headers
  766.     (a, b) Bigint *a, *b;
  767. #else
  768.     (Bigint *a, Bigint *b)
  769. #endif
  770. {
  771.     unsigned Long *xa, *xa0, *xb, *xb0;
  772.     int i, j;
  773.  
  774.     i = a->wds;
  775.     j = b->wds;
  776. #ifdef DEBUG
  777.     if (i > 1 && !a->x[i-1])
  778.         Bug("cmp called with a->x[a->wds-1] == 0");
  779.     if (j > 1 && !b->x[j-1])
  780.         Bug("cmp called with b->x[b->wds-1] == 0");
  781. #endif
  782.     if (i -= j)
  783.         return i;
  784.     xa0 = a->x;
  785.     xa = xa0 + j;
  786.     xb0 = b->x;
  787.     xb = xb0 + j;
  788.     for(;;) {
  789.         if (*--xa != *--xb)
  790.             return *xa < *xb ? -1 : 1;
  791.         if (xa <= xa0)
  792.             break;
  793.         }
  794.     return 0;
  795.     }
  796.  
  797.  static Bigint *
  798. diff
  799. #ifdef KR_headers
  800.     (a, b) Bigint *a, *b;
  801. #else
  802.     (Bigint *a, Bigint *b)
  803. #endif
  804. {
  805.     Bigint *c;
  806.     int i, wa, wb;
  807.     Long borrow, y;    /* We need signed shifts here. */
  808.     unsigned Long *xa, *xae, *xb, *xbe, *xc;
  809. #ifdef Pack_32
  810.     Long z;
  811. #endif
  812.  
  813.     i = cmp(a,b);
  814.     if (!i) {
  815.         c = Balloc(0);
  816.         c->wds = 1;
  817.         c->x[0] = 0;
  818.         return c;
  819.         }
  820.     if (i < 0) {
  821.         c = a;
  822.         a = b;
  823.         b = c;
  824.         i = 1;
  825.         }
  826.     else
  827.         i = 0;
  828.     c = Balloc(a->k);
  829.     c->sign = i;
  830.     wa = a->wds;
  831.     xa = a->x;
  832.     xae = xa + wa;
  833.     wb = b->wds;
  834.     xb = b->x;
  835.     xbe = xb + wb;
  836.     xc = c->x;
  837.     borrow = 0;
  838. #ifdef Pack_32
  839.     do {
  840.         y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
  841.         borrow = y >> 16;
  842.         Sign_Extend(borrow, y);
  843.         z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
  844.         borrow = z >> 16;
  845.         Sign_Extend(borrow, z);
  846.         Storeinc(xc, z, y);
  847.         }
  848.         while(xb < xbe);
  849.     while(xa < xae) {
  850.         y = (*xa & 0xffff) + borrow;
  851.         borrow = y >> 16;
  852.         Sign_Extend(borrow, y);
  853.         z = (*xa++ >> 16) + borrow;
  854.         borrow = z >> 16;
  855.         Sign_Extend(borrow, z);
  856.         Storeinc(xc, z, y);
  857.         }
  858. #else
  859.     do {
  860.         y = *xa++ - *xb++ + borrow;
  861.         borrow = y >> 16;
  862.         Sign_Extend(borrow, y);
  863.         *xc++ = y & 0xffff;
  864.         }
  865.         while(xb < xbe);
  866.     while(xa < xae) {
  867.         y = *xa++ + borrow;
  868.         borrow = y >> 16;
  869.         Sign_Extend(borrow, y);
  870.         *xc++ = y & 0xffff;
  871.         }
  872. #endif
  873.     while(!*--xc)
  874.         wa--;
  875.     c->wds = wa;
  876.     return c;
  877.     }
  878.  
  879.  static double
  880. ulp
  881. #ifdef KR_headers
  882.     (x) double x;
  883. #else
  884.     (double x)
  885. #endif
  886. {
  887.     register Long L;
  888.     double a;
  889.  
  890.     L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
  891. #ifndef Sudden_Underflow
  892.     if (L > 0) {
  893. #endif
  894. #ifdef IBM
  895.         L |= Exp_msk1 >> 4;
  896. #endif
  897.         word0(a) = L;
  898.         word1(a) = 0;
  899. #ifndef Sudden_Underflow
  900.         }
  901.     else {
  902.         L = -L >> Exp_shift;
  903.         if (L < Exp_shift) {
  904.             word0(a) = 0x80000 >> L;
  905.             word1(a) = 0;
  906.             }
  907.         else {
  908.             word0(a) = 0;
  909.             L -= Exp_shift;
  910.             word1(a) = L >= 31 ? 1 : 1 << 31 - L;
  911.             }
  912.         }
  913. #endif
  914.     return a;
  915.     }
  916.  
  917.  static double
  918. b2d
  919. #ifdef KR_headers
  920.     (a, e) Bigint *a; int *e;
  921. #else
  922.     (Bigint *a, int *e)
  923. #endif
  924. {
  925.     unsigned Long *xa, *xa0, w, y, z;
  926.     int k;
  927.     double d;
  928. #ifdef VAX
  929.     unsigned Long d0, d1;
  930. #else
  931. #define d0 word0(d)
  932. #define d1 word1(d)
  933. #endif
  934.  
  935.     xa0 = a->x;
  936.     xa = xa0 + a->wds;
  937.     y = *--xa;
  938. #ifdef DEBUG
  939.     if (!y) Bug("zero y in b2d");
  940. #endif
  941.     k = hi0bits(y);
  942.     *e = 32 - k;
  943. #ifdef Pack_32
  944.     if (k < Ebits) {
  945.         d0 = Exp_1 | y >> Ebits - k;
  946.         w = xa > xa0 ? *--xa : 0;
  947.         d1 = y << (32-Ebits) + k | w >> Ebits - k;
  948.         goto ret_d;
  949.         }
  950.     z = xa > xa0 ? *--xa : 0;
  951.     if (k -= Ebits) {
  952.         d0 = Exp_1 | y << k | z >> 32 - k;
  953.         y = xa > xa0 ? *--xa : 0;
  954.         d1 = z << k | y >> 32 - k;
  955.         }
  956.     else {
  957.         d0 = Exp_1 | y;
  958.         d1 = z;
  959.         }
  960. #else
  961.     if (k < Ebits + 16) {
  962.         z = xa > xa0 ? *--xa : 0;
  963.         d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
  964.         w = xa > xa0 ? *--xa : 0;
  965.         y = xa > xa0 ? *--xa : 0;
  966.         d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
  967.         goto ret_d;
  968.         }
  969.     z = xa > xa0 ? *--xa : 0;
  970.     w = xa > xa0 ? *--xa : 0;
  971.     k -= Ebits + 16;
  972.     d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
  973.     y = xa > xa0 ? *--xa : 0;
  974.     d1 = w << k + 16 | y << k;
  975. #endif
  976.  ret_d:
  977. #ifdef VAX
  978.     word0(d) = d0 >> 16 | d0 << 16;
  979.     word1(d) = d1 >> 16 | d1 << 16;
  980. #else
  981. #undef d0
  982. #undef d1
  983. #endif
  984.     return d;
  985.     }
  986.  
  987.  static Bigint *
  988. d2b
  989. #ifdef KR_headers
  990.     (d, e, bits) double d; int *e, *bits;
  991. #else
  992.     (double d, int *e, int *bits)
  993. #endif
  994. {
  995.     Bigint *b;
  996.     int de, i, k;
  997.     unsigned Long *x, y, z;
  998. #ifdef VAX
  999.     unsigned Long d0, d1;
  1000.     d0 = word0(d) >> 16 | word0(d) << 16;
  1001.     d1 = word1(d) >> 16 | word1(d) << 16;
  1002. #else
  1003. #define d0 word0(d)
  1004. #define d1 word1(d)
  1005. #endif
  1006.  
  1007. #ifdef Pack_32
  1008.     b = Balloc(1);
  1009. #else
  1010.     b = Balloc(2);
  1011. #endif
  1012.     x = b->x;
  1013.  
  1014.     z = d0 & Frac_mask;
  1015.     d0 &= 0x7fffffff;    /* clear sign bit, which we ignore */
  1016. #ifdef Sudden_Underflow
  1017.     de = (int)(d0 >> Exp_shift);
  1018. #ifndef IBM
  1019.     z |= Exp_msk11;
  1020. #endif
  1021. #else
  1022.     if (de = (int)(d0 >> Exp_shift))
  1023.         z |= Exp_msk1;
  1024. #endif
  1025. #ifdef Pack_32
  1026.     if (y = d1) {
  1027.         if (k = lo0bits(&y)) {
  1028.             x[0] = y | z << 32 - k;
  1029.             z >>= k;
  1030.             }
  1031.         else
  1032.             x[0] = y;
  1033.         i = b->wds = (x[1] = z) ? 2 : 1;
  1034.         }
  1035.     else {
  1036. #ifdef DEBUG
  1037.         if (!z)
  1038.             Bug("Zero passed to d2b");
  1039. #endif
  1040.         k = lo0bits(&z);
  1041.         x[0] = z;
  1042.         i = b->wds = 1;
  1043.         k += 32;
  1044.         }
  1045. #else
  1046.     if (y = d1) {
  1047.         if (k = lo0bits(&y))
  1048.             if (k >= 16) {
  1049.                 x[0] = y | z << 32 - k & 0xffff;
  1050.                 x[1] = z >> k - 16 & 0xffff;
  1051.                 x[2] = z >> k;
  1052.                 i = 2;
  1053.                 }
  1054.             else {
  1055.                 x[0] = y & 0xffff;
  1056.                 x[1] = y >> 16 | z << 16 - k & 0xffff;
  1057.                 x[2] = z >> k & 0xffff;
  1058.                 x[3] = z >> k+16;
  1059.                 i = 3;
  1060.                 }
  1061.         else {
  1062.             x[0] = y & 0xffff;
  1063.             x[1] = y >> 16;
  1064.             x[2] = z & 0xffff;
  1065.             x[3] = z >> 16;
  1066.             i = 3;
  1067.             }
  1068.         }
  1069.     else {
  1070. #ifdef DEBUG
  1071.         if (!z)
  1072.             Bug("Zero passed to d2b");
  1073. #endif
  1074.         k = lo0bits(&z);
  1075.         if (k >= 16) {
  1076.             x[0] = z;
  1077.             i = 0;
  1078.             }
  1079.         else {
  1080.             x[0] = z & 0xffff;
  1081.             x[1] = z >> 16;
  1082.             i = 1;
  1083.             }
  1084.         k += 32;
  1085.         }
  1086.     while(!x[i])
  1087.         --i;
  1088.     b->wds = i + 1;
  1089. #endif
  1090. #ifndef Sudden_Underflow
  1091.     if (de) {
  1092. #endif
  1093. #ifdef IBM
  1094.         *e = (de - Bias - (P-1) << 2) + k;
  1095.         *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
  1096. #else
  1097.         *e = de - Bias - (P-1) + k;
  1098.         *bits = P - k;
  1099. #endif
  1100. #ifndef Sudden_Underflow
  1101.         }
  1102.     else {
  1103.         *e = de - Bias - (P-1) + 1 + k;
  1104. #ifdef Pack_32
  1105.         *bits = 32*i - hi0bits(x[i-1]);
  1106. #else
  1107.         *bits = (i+2)*16 - hi0bits(x[i]);
  1108. #endif
  1109.         }
  1110. #endif
  1111.     return b;
  1112.     }
  1113. #undef d0
  1114. #undef d1
  1115.  
  1116.  static double
  1117. ratio
  1118. #ifdef KR_headers
  1119.     (a, b) Bigint *a, *b;
  1120. #else
  1121.     (Bigint *a, Bigint *b)
  1122. #endif
  1123. {
  1124.     double da, db;
  1125.     int k, ka, kb;
  1126.  
  1127.     da = b2d(a, &ka);
  1128.     db = b2d(b, &kb);
  1129. #ifdef Pack_32
  1130.     k = ka - kb + 32*(a->wds - b->wds);
  1131. #else
  1132.     k = ka - kb + 16*(a->wds - b->wds);
  1133. #endif
  1134. #ifdef IBM
  1135.     if (k > 0) {
  1136.         word0(da) += (k >> 2)*Exp_msk1;
  1137.         if (k &= 3)
  1138.             da *= 1 << k;
  1139.         }
  1140.     else {
  1141.         k = -k;
  1142.         word0(db) += (k >> 2)*Exp_msk1;
  1143.         if (k &= 3)
  1144.             db *= 1 << k;
  1145.         }
  1146. #else
  1147.     if (k > 0)
  1148.         word0(da) += k*Exp_msk1;
  1149.     else {
  1150.         k = -k;
  1151.         word0(db) += k*Exp_msk1;
  1152.         }
  1153. #endif
  1154.     return da / db;
  1155.     }
  1156.  
  1157.  static double
  1158. tens[] = {
  1159.         1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
  1160.         1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
  1161.         1e20, 1e21, 1e22
  1162. #ifdef VAX
  1163.         , 1e23, 1e24
  1164. #endif
  1165.         };
  1166.  
  1167.  static double
  1168. #ifdef IEEE_Arith
  1169. bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
  1170. static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
  1171. #define n_bigtens 5
  1172. #else
  1173. #ifdef IBM
  1174. bigtens[] = { 1e16, 1e32, 1e64 };
  1175. static double tinytens[] = { 1e-16, 1e-32, 1e-64 };
  1176. #define n_bigtens 3
  1177. #else
  1178. bigtens[] = { 1e16, 1e32 };
  1179. static double tinytens[] = { 1e-16, 1e-32 };
  1180. #define n_bigtens 2
  1181. #endif
  1182. #endif
  1183.  
  1184.  double
  1185. strtod
  1186. #ifdef KR_headers
  1187.     (s00, se) CONST char *s00; char **se;
  1188. #else
  1189.     (CONST char *s00, char **se)
  1190. #endif
  1191. {
  1192.     int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
  1193.          e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
  1194.     CONST char *s, *s0, *s1;
  1195.     double aadj, aadj1, adj, rv, rv0;
  1196.     Long L;
  1197.     unsigned Long y, z;
  1198.     Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
  1199.     sign = nz0 = nz = 0;
  1200.     rv = 0.;
  1201.     for(s = s00;;s++) switch(*s) {
  1202.         case '-':
  1203.             sign = 1;
  1204.             /* no break */
  1205.         case '+':
  1206.             if (*++s)
  1207.                 goto break2;
  1208.             /* no break */
  1209.         case 0:
  1210.             s = s00;
  1211.             goto ret;
  1212.         case '\t':
  1213.         case '\n':
  1214.         case '\v':
  1215.         case '\f':
  1216.         case '\r':
  1217.         case ' ':
  1218.             continue;
  1219.         default:
  1220.             goto break2;
  1221.         }
  1222.  break2:
  1223.     if (*s == '0') {
  1224.         nz0 = 1;
  1225.         while(*++s == '0') ;
  1226.         if (!*s)
  1227.             goto ret;
  1228.         }
  1229.     s0 = s;
  1230.     y = z = 0;
  1231.     for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
  1232.         if (nd < 9)
  1233.             y = 10*y + c - '0';
  1234.         else if (nd < 16)
  1235.             z = 10*z + c - '0';
  1236.     nd0 = nd;
  1237.     if (c == '.') {
  1238.         c = *++s;
  1239.         if (!nd) {
  1240.             for(; c == '0'; c = *++s)
  1241.                 nz++;
  1242.             if (c > '0' && c <= '9') {
  1243.                 s0 = s;
  1244.                 nf += nz;
  1245.                 nz = 0;
  1246.                 goto have_dig;
  1247.                 }
  1248.             goto dig_done;
  1249.             }
  1250.         for(; c >= '0' && c <= '9'; c = *++s) {
  1251.  have_dig:
  1252.             nz++;
  1253.             if (c -= '0') {
  1254.                 nf += nz;
  1255.                 for(i = 1; i < nz; i++)
  1256.                     if (nd++ < 9)
  1257.                         y *= 10;
  1258.                     else if (nd <= DBL_DIG + 1)
  1259.                         z *= 10;
  1260.                 if (nd++ < 9)
  1261.                     y = 10*y + c;
  1262.                 else if (nd <= DBL_DIG + 1)
  1263.                     z = 10*z + c;
  1264.                 nz = 0;
  1265.                 }
  1266.             }
  1267.         }
  1268.  dig_done:
  1269.     e = 0;
  1270.     if (c == 'e' || c == 'E') {
  1271.         if (!nd && !nz && !nz0) {
  1272.             s = s00;
  1273.             goto ret;
  1274.             }
  1275.         s00 = s;
  1276.         esign = 0;
  1277.         switch(c = *++s) {
  1278.             case '-':
  1279.                 esign = 1;
  1280.             case '+':
  1281.                 c = *++s;
  1282.             }
  1283.         if (c >= '0' && c <= '9') {
  1284.             while(c == '0')
  1285.                 c = *++s;
  1286.             if (c > '0' && c <= '9') {
  1287.                 L = c - '0';
  1288.                 s1 = s;
  1289.                 while((c = *++s) >= '0' && c <= '9')
  1290.                     L = 10*L + c - '0';
  1291.                 if (s - s1 > 8 || L > 19999)
  1292.                     /* Avoid confusion from exponents
  1293.                      * so large that e might overflow.
  1294.                      */
  1295.                     e = 19999; /* safe for 16 bit ints */
  1296.                 else
  1297.                     e = (int)L;
  1298.                 if (esign)
  1299.                     e = -e;
  1300.                 }
  1301.             else
  1302.                 e = 0;
  1303.             }
  1304.         else
  1305.             s = s00;
  1306.         }
  1307.     if (!nd) {
  1308.         if (!nz && !nz0)
  1309.             s = s00;
  1310.         goto ret;
  1311.         }
  1312.     e1 = e -= nf;
  1313.  
  1314.     /* Now we have nd0 digits, starting at s0, followed by a
  1315.      * decimal point, followed by nd-nd0 digits.  The number we're
  1316.      * after is the integer represented by those digits times
  1317.      * 10**e */
  1318.  
  1319.     if (!nd0)
  1320.         nd0 = nd;
  1321.     k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
  1322.     rv = y;
  1323.     if (k > 9)
  1324.         rv = tens[k - 9] * rv + z;
  1325.     bd0 = 0;
  1326.     if (nd <= DBL_DIG
  1327. #ifndef RND_PRODQUOT
  1328.         && FLT_ROUNDS == 1
  1329. #endif
  1330.             ) {
  1331.         if (!e)
  1332.             goto ret;
  1333.         if (e > 0) {
  1334.             if (e <= Ten_pmax) {
  1335. #ifdef VAX
  1336.                 goto vax_ovfl_check;
  1337. #else
  1338.                 /* rv = */ rounded_product(rv, tens[e]);
  1339.                 goto ret;
  1340. #endif
  1341.                 }
  1342.             i = DBL_DIG - nd;
  1343.             if (e <= Ten_pmax + i) {
  1344.                 /* A fancier test would sometimes let us do
  1345.                  * this for larger i values.
  1346.                  */
  1347.                 e -= i;
  1348.                 rv *= tens[i];
  1349. #ifdef VAX
  1350.                 /* VAX exponent range is so narrow we must
  1351.                  * worry about overflow here...
  1352.                  */
  1353.  vax_ovfl_check:
  1354.                 word0(rv) -= P*Exp_msk1;
  1355.                 /* rv = */ rounded_product(rv, tens[e]);
  1356.                 if ((word0(rv) & Exp_mask)
  1357.                  > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
  1358.                     goto ovfl;
  1359.                 word0(rv) += P*Exp_msk1;
  1360. #else
  1361.                 /* rv = */ rounded_product(rv, tens[e]);
  1362. #endif
  1363.                 goto ret;
  1364.                 }
  1365.             }
  1366. #ifndef Inaccurate_Divide
  1367.         else if (e >= -Ten_pmax) {
  1368.             /* rv = */ rounded_quotient(rv, tens[-e]);
  1369.             goto ret;
  1370.             }
  1371. #endif
  1372.         }
  1373.     e1 += nd - k;
  1374.  
  1375.     /* Get starting approximation = rv * 10**e1 */
  1376.  
  1377.     if (e1 > 0) {
  1378.         if (i = e1 & 15)
  1379.             rv *= tens[i];
  1380.         if (e1 &= ~15) {
  1381.             if (e1 > DBL_MAX_10_EXP) {
  1382.  ovfl:
  1383.                 errno = ERANGE;
  1384. #ifdef __STDC__
  1385.                 rv = HUGE_VAL;
  1386. #else
  1387.                 /* Can't trust HUGE_VAL */
  1388. #ifdef IEEE_Arith
  1389.                 word0(rv) = Exp_mask;
  1390.                 word1(rv) = 0;
  1391. #else
  1392.                 word0(rv) = Big0;
  1393.                 word1(rv) = Big1;
  1394. #endif
  1395. #endif
  1396.                 if (bd0)
  1397.                     goto retfree;
  1398.                 goto ret;
  1399.                 }
  1400.             if (e1 >>= 4) {
  1401.                 for(j = 0; e1 > 1; j++, e1 >>= 1)
  1402.                     if (e1 & 1)
  1403.                         rv *= bigtens[j];
  1404.             /* The last multiplication could overflow. */
  1405.                 word0(rv) -= P*Exp_msk1;
  1406.                 rv *= bigtens[j];
  1407.                 if ((z = word0(rv) & Exp_mask)
  1408.                  > Exp_msk1*(DBL_MAX_EXP+Bias-P))
  1409.                     goto ovfl;
  1410.                 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
  1411.                     /* set to largest number */
  1412.                     /* (Can't trust DBL_MAX) */
  1413.                     word0(rv) = Big0;
  1414.                     word1(rv) = Big1;
  1415.                     }
  1416.                 else
  1417.                     word0(rv) += P*Exp_msk1;
  1418.                 }
  1419.  
  1420.             }
  1421.         }
  1422.     else if (e1 < 0) {
  1423.         e1 = -e1;
  1424.         if (i = e1 & 15)
  1425.             rv /= tens[i];
  1426.         if (e1 &= ~15) {
  1427.             e1 >>= 4;
  1428.             if (e1 >= 1 << n_bigtens)
  1429.                 goto undfl;
  1430.             for(j = 0; e1 > 1; j++, e1 >>= 1)
  1431.                 if (e1 & 1)
  1432.                     rv *= tinytens[j];
  1433.             /* The last multiplication could underflow. */
  1434.             rv0 = rv;
  1435.             rv *= tinytens[j];
  1436.             if (!rv) {
  1437.                 rv = 2.*rv0;
  1438.                 rv *= tinytens[j];
  1439.                 if (!rv) {
  1440.  undfl:
  1441.                     rv = 0.;
  1442.                     errno = ERANGE;
  1443.                     if (bd0)
  1444.                         goto retfree;
  1445.                     goto ret;
  1446.                     }
  1447.                 word0(rv) = Tiny0;
  1448.                 word1(rv) = Tiny1;
  1449.                 /* The refinement below will clean
  1450.                  * this approximation up.
  1451.                  */
  1452.                 }
  1453.             }
  1454.         }
  1455.  
  1456.     /* Now the hard part -- adjusting rv to the correct value.*/
  1457.  
  1458.     /* Put digits into bd: true value = bd * 10^e */
  1459.  
  1460.     bd0 = s2b(s0, nd0, nd, y);
  1461.  
  1462.     for(;;) {
  1463.         bd = Balloc(bd0->k);
  1464.         Bcopy(bd, bd0);
  1465.         bb = d2b(rv, &bbe, &bbbits);    /* rv = bb * 2^bbe */
  1466.         bs = i2b(1);
  1467.  
  1468.         if (e >= 0) {
  1469.             bb2 = bb5 = 0;
  1470.             bd2 = bd5 = e;
  1471.             }
  1472.         else {
  1473.             bb2 = bb5 = -e;
  1474.             bd2 = bd5 = 0;
  1475.             }
  1476.         if (bbe >= 0)
  1477.             bb2 += bbe;
  1478.         else
  1479.             bd2 -= bbe;
  1480.         bs2 = bb2;
  1481. #ifdef Sudden_Underflow
  1482. #ifdef IBM
  1483.         j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
  1484. #else
  1485.         j = P + 1 - bbbits;
  1486. #endif
  1487. #else
  1488.         i = bbe + bbbits - 1;    /* logb(rv) */
  1489.         if (i < Emin)    /* denormal */
  1490.             j = bbe + (P-Emin);
  1491.         else
  1492.             j = P + 1 - bbbits;
  1493. #endif
  1494.         bb2 += j;
  1495.         bd2 += j;
  1496.         i = bb2 < bd2 ? bb2 : bd2;
  1497.         if (i > bs2)
  1498.             i = bs2;
  1499.         if (i > 0) {
  1500.             bb2 -= i;
  1501.             bd2 -= i;
  1502.             bs2 -= i;
  1503.             }
  1504.         if (bb5 > 0) {
  1505.             bs = pow5mult(bs, bb5);
  1506.             bb1 = mult(bs, bb);
  1507.             Bfree(bb);
  1508.             bb = bb1;
  1509.             }
  1510.         if (bb2 > 0)
  1511.             bb = lshift(bb, bb2);
  1512.         if (bd5 > 0)
  1513.             bd = pow5mult(bd, bd5);
  1514.         if (bd2 > 0)
  1515.             bd = lshift(bd, bd2);
  1516.         if (bs2 > 0)
  1517.             bs = lshift(bs, bs2);
  1518.         delta = diff(bb, bd);
  1519.         dsign = delta->sign;
  1520.         delta->sign = 0;
  1521.         i = cmp(delta, bs);
  1522.         if (i < 0) {
  1523.             /* Error is less than half an ulp -- check for
  1524.              * special case of mantissa a power of two.
  1525.              */
  1526.             if (dsign || word1(rv) || word0(rv) & Bndry_mask)
  1527.                 break;
  1528.             delta = lshift(delta,Log2P);
  1529.             if (cmp(delta, bs) > 0)
  1530.                 goto drop_down;
  1531.             break;
  1532.             }
  1533.         if (i == 0) {
  1534.             /* exactly half-way between */
  1535.             if (dsign) {
  1536.                 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
  1537.                  &&  word1(rv) == 0xffffffff) {
  1538.                     /*boundary case -- increment exponent*/
  1539.                     word0(rv) = (word0(rv) & Exp_mask)
  1540.                         + Exp_msk1
  1541. #ifdef IBM
  1542.                         | Exp_msk1 >> 4
  1543. #endif
  1544.                         ;
  1545.                     word1(rv) = 0;
  1546.                     break;
  1547.                     }
  1548.                 }
  1549.             else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
  1550.  drop_down:
  1551.                 /* boundary case -- decrement exponent */
  1552. #ifdef Sudden_Underflow
  1553.                 L = word0(rv) & Exp_mask;
  1554. #ifdef IBM
  1555.                 if (L <  Exp_msk1)
  1556. #else
  1557.                 if (L <= Exp_msk1)
  1558. #endif
  1559.                     goto undfl;
  1560.                 L -= Exp_msk1;
  1561. #else
  1562.                 L = (word0(rv) & Exp_mask) - Exp_msk1;
  1563. #endif
  1564.                 word0(rv) = L | Bndry_mask1;
  1565.                 word1(rv) = 0xffffffff;
  1566. #ifdef IBM
  1567.                 goto cont;
  1568. #else
  1569.                 break;
  1570. #endif
  1571.                 }
  1572. #ifndef ROUND_BIASED
  1573.             if (!(word1(rv) & LSB))
  1574.                 break;
  1575. #endif
  1576.             if (dsign)
  1577.                 rv += ulp(rv);
  1578. #ifndef ROUND_BIASED
  1579.             else {
  1580.                 rv -= ulp(rv);
  1581. #ifndef Sudden_Underflow
  1582.                 if (!rv)
  1583.                     goto undfl;
  1584. #endif
  1585.                 }
  1586. #endif
  1587.             break;
  1588.             }
  1589.         if ((aadj = ratio(delta, bs)) <= 2.) {
  1590.             if (dsign)
  1591.                 aadj = aadj1 = 1.;
  1592.             else if (word1(rv) || word0(rv) & Bndry_mask) {
  1593. #ifndef Sudden_Underflow
  1594.                 if (word1(rv) == Tiny1 && !word0(rv))
  1595.                     goto undfl;
  1596. #endif
  1597.                 aadj = 1.;
  1598.                 aadj1 = -1.;
  1599.                 }
  1600.             else {
  1601.                 /* special case -- power of FLT_RADIX to be */
  1602.                 /* rounded down... */
  1603.  
  1604.                 if (aadj < 2./FLT_RADIX)
  1605.                     aadj = 1./FLT_RADIX;
  1606.                 else
  1607.                     aadj *= 0.5;
  1608.                 aadj1 = -aadj;
  1609.                 }
  1610.             }
  1611.         else {
  1612.             aadj *= 0.5;
  1613.             aadj1 = dsign ? aadj : -aadj;
  1614. #ifdef Check_FLT_ROUNDS
  1615.             switch(FLT_ROUNDS) {
  1616.                 case 2: /* towards +infinity */
  1617.                     aadj1 -= 0.5;
  1618.                     break;
  1619.                 case 0: /* towards 0 */
  1620.                 case 3: /* towards -infinity */
  1621.                     aadj1 += 0.5;
  1622.                 }
  1623. #else
  1624.             if (FLT_ROUNDS == 0)
  1625.                 aadj1 += 0.5;
  1626. #endif
  1627.             }
  1628.         y = word0(rv) & Exp_mask;
  1629.  
  1630.         /* Check for overflow */
  1631.  
  1632.         if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
  1633.             rv0 = rv;
  1634.             word0(rv) -= P*Exp_msk1;
  1635.             adj = aadj1 * ulp(rv);
  1636.             rv += adj;
  1637.             if ((word0(rv) & Exp_mask) >=
  1638.                     Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
  1639.                 if (word0(rv0) == Big0 && word1(rv0) == Big1)
  1640.                     goto ovfl;
  1641.                 word0(rv) = Big0;
  1642.                 word1(rv) = Big1;
  1643.                 goto cont;
  1644.                 }
  1645.             else
  1646.                 word0(rv) += P*Exp_msk1;
  1647.             }
  1648.         else {
  1649. #ifdef Sudden_Underflow
  1650.             if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
  1651.                 rv0 = rv;
  1652.                 word0(rv) += P*Exp_msk1;
  1653.                 adj = aadj1 * ulp(rv);
  1654.                 rv += adj;
  1655. #ifdef IBM
  1656.                 if ((word0(rv) & Exp_mask) <  P*Exp_msk1)
  1657. #else
  1658.                 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
  1659. #endif
  1660.                     {
  1661.                     if (word0(rv0) == Tiny0
  1662.                      && word1(rv0) == Tiny1)
  1663.                         goto undfl;
  1664.                     word0(rv) = Tiny0;
  1665.                     word1(rv) = Tiny1;
  1666.                     goto cont;
  1667.                     }
  1668.                 else
  1669.                     word0(rv) -= P*Exp_msk1;
  1670.                 }
  1671.             else {
  1672.                 adj = aadj1 * ulp(rv);
  1673.                 rv += adj;
  1674.                 }
  1675. #else
  1676.             /* Compute adj so that the IEEE rounding rules will
  1677.              * correctly round rv + adj in some half-way cases.
  1678.              * If rv * ulp(rv) is denormalized (i.e.,
  1679.              * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
  1680.              * trouble from bits lost to denormalization;
  1681.              * example: 1.2e-307 .
  1682.              */
  1683.             if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
  1684.                 aadj1 = (double)(int)(aadj + 0.5);
  1685.                 if (!dsign)
  1686.                     aadj1 = -aadj1;
  1687.                 }
  1688.             adj = aadj1 * ulp(rv);
  1689.             rv += adj;
  1690. #endif
  1691.             }
  1692.         z = word0(rv) & Exp_mask;
  1693.         if (y == z) {
  1694.             /* Can we stop now? */
  1695.             L = aadj;
  1696.             aadj -= L;
  1697.             /* The tolerances below are conservative. */
  1698.             if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
  1699.                 if (aadj < .4999999 || aadj > .5000001)
  1700.                     break;
  1701.                 }
  1702.             else if (aadj < .4999999/FLT_RADIX)
  1703.                 break;
  1704.             }
  1705.  cont:
  1706.         Bfree(bb);
  1707.         Bfree(bd);
  1708.         Bfree(bs);
  1709.         Bfree(delta);
  1710.         }
  1711.  retfree:
  1712.     Bfree(bb);
  1713.     Bfree(bd);
  1714.     Bfree(bs);
  1715.     Bfree(bd0);
  1716.     Bfree(delta);
  1717.  ret:
  1718.     if (se)
  1719.         *se = (char *)s;
  1720.     return sign ? -rv : rv;
  1721.     }
  1722.  
  1723.  static int
  1724. quorem
  1725. #ifdef KR_headers
  1726.     (b, S) Bigint *b, *S;
  1727. #else
  1728.     (Bigint *b, Bigint *S)
  1729. #endif
  1730. {
  1731.     int n;
  1732.     Long borrow, y;
  1733.     unsigned Long carry, q, ys;
  1734.     unsigned Long *bx, *bxe, *sx, *sxe;
  1735. #ifdef Pack_32
  1736.     Long z;
  1737.     unsigned Long si, zs;
  1738. #endif
  1739.  
  1740.     n = S->wds;
  1741. #ifdef DEBUG
  1742.     /*debug*/ if (b->wds > n)
  1743.     /*debug*/    Bug("oversize b in quorem");
  1744. #endif
  1745.     if (b->wds < n)
  1746.         return 0;
  1747.     sx = S->x;
  1748.     sxe = sx + --n;
  1749.     bx = b->x;
  1750.     bxe = bx + n;
  1751.     q = *bxe / (*sxe + 1);    /* ensure q <= true quotient */
  1752. #ifdef DEBUG
  1753.     /*debug*/ if (q > 9)
  1754.     /*debug*/    Bug("oversized quotient in quorem");
  1755. #endif
  1756.     if (q) {
  1757.         borrow = 0;
  1758.         carry = 0;
  1759.         do {
  1760. #ifdef Pack_32
  1761.             si = *sx++;
  1762.             ys = (si & 0xffff) * q + carry;
  1763.             zs = (si >> 16) * q + (ys >> 16);
  1764.             carry = zs >> 16;
  1765.             y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
  1766.             borrow = y >> 16;
  1767.             Sign_Extend(borrow, y);
  1768.             z = (*bx >> 16) - (zs & 0xffff) + borrow;
  1769.             borrow = z >> 16;
  1770.             Sign_Extend(borrow, z);
  1771.             Storeinc(bx, z, y);
  1772. #else
  1773.             ys = *sx++ * q + carry;
  1774.             carry = ys >> 16;
  1775.             y = *bx - (ys & 0xffff) + borrow;
  1776.             borrow = y >> 16;
  1777.             Sign_Extend(borrow, y);
  1778.             *bx++ = y & 0xffff;
  1779. #endif
  1780.             }
  1781.             while(sx <= sxe);
  1782.         if (!*bxe) {
  1783.             bx = b->x;
  1784.             while(--bxe > bx && !*bxe)
  1785.                 --n;
  1786.             b->wds = n;
  1787.             }
  1788.         }
  1789.     if (cmp(b, S) >= 0) {
  1790.         q++;
  1791.         borrow = 0;
  1792.         carry = 0;
  1793.         bx = b->x;
  1794.         sx = S->x;
  1795.         do {
  1796. #ifdef Pack_32
  1797.             si = *sx++;
  1798.             ys = (si & 0xffff) + carry;
  1799.             zs = (si >> 16) + (ys >> 16);
  1800.             carry = zs >> 16;
  1801.             y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
  1802.             borrow = y >> 16;
  1803.             Sign_Extend(borrow, y);
  1804.             z = (*bx >> 16) - (zs & 0xffff) + borrow;
  1805.             borrow = z >> 16;
  1806.             Sign_Extend(borrow, z);
  1807.             Storeinc(bx, z, y);
  1808. #else
  1809.             ys = *sx++ + carry;
  1810.             carry = ys >> 16;
  1811.             y = *bx - (ys & 0xffff) + borrow;
  1812.             borrow = y >> 16;
  1813.             Sign_Extend(borrow, y);
  1814.             *bx++ = y & 0xffff;
  1815. #endif
  1816.             }
  1817.             while(sx <= sxe);
  1818.         bx = b->x;
  1819.         bxe = bx + n;
  1820.         if (!*bxe) {
  1821.             while(--bxe > bx && !*bxe)
  1822.                 --n;
  1823.             b->wds = n;
  1824.             }
  1825.         }
  1826.     return q;
  1827.     }
  1828.  
  1829. /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
  1830.  *
  1831.  * Inspired by "How to Print Floating-Point Numbers Accurately" by
  1832.  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
  1833.  *
  1834.  * Modifications:
  1835.  *    1. Rather than iterating, we use a simple numeric overestimate
  1836.  *       to determine k = floor(log10(d)).  We scale relevant
  1837.  *       quantities using O(log2(k)) rather than O(k) multiplications.
  1838.  *    2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
  1839.  *       try to generate digits strictly left to right.  Instead, we
  1840.  *       compute with fewer bits and propagate the carry if necessary
  1841.  *       when rounding the final digit up.  This is often faster.
  1842.  *    3. Under the assumption that input will be rounded nearest,
  1843.  *       mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
  1844.  *       That is, we allow equality in stopping tests when the
  1845.  *       round-nearest rule will give the same floating-point value
  1846.  *       as would satisfaction of the stopping test with strict
  1847.  *       inequality.
  1848.  *    4. We remove common factors of powers of 2 from relevant
  1849.  *       quantities.
  1850.  *    5. When converting floating-point integers less than 1e16,
  1851.  *       we use floating-point arithmetic rather than resorting
  1852.  *       to multiple-precision integers.
  1853.  *    6. When asked to produce fewer than 15 digits, we first try
  1854.  *       to get by with floating-point arithmetic; we resort to
  1855.  *       multiple-precision integer arithmetic only if we cannot
  1856.  *       guarantee that the floating-point calculation has given
  1857.  *       the correctly rounded result.  For k requested digits and
  1858.  *       "uniformly" distributed input, the probability is
  1859.  *       something like 10^(k-15) that we must resort to the Long
  1860.  *       calculation.
  1861.  */
  1862.  
  1863.  char *
  1864. dtoa
  1865. #ifdef KR_headers
  1866.     (d, mode, ndigits, decpt, sign, rve)
  1867.     double d; int mode, ndigits, *decpt, *sign; char **rve;
  1868. #else
  1869.     (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
  1870. #endif
  1871. {
  1872.  /*    Arguments ndigits, decpt, sign are similar to those
  1873.     of ecvt and fcvt; trailing zeros are suppressed from
  1874.     the returned string.  If not null, *rve is set to point
  1875.     to the end of the return value.  If d is +-Infinity or NaN,
  1876.     then *decpt is set to 9999.
  1877.  
  1878.     mode:
  1879.         0 ==> shortest string that yields d when read in
  1880.             and rounded to nearest.
  1881.         1 ==> like 0, but with Steele & White stopping rule;
  1882.             e.g. with IEEE P754 arithmetic , mode 0 gives
  1883.             1e23 whereas mode 1 gives 9.999999999999999e22.
  1884.         2 ==> max(1,ndigits) significant digits.  This gives a
  1885.             return value similar to that of ecvt, except
  1886.             that trailing zeros are suppressed.
  1887.         3 ==> through ndigits past the decimal point.  This
  1888.             gives a return value similar to that from fcvt,
  1889.             except that trailing zeros are suppressed, and
  1890.             ndigits can be negative.
  1891.         4-9 should give the same return values as 2-3, i.e.,
  1892.             4 <= mode <= 9 ==> same return as mode
  1893.             2 + (mode & 1).  These modes are mainly for
  1894.             debugging; often they run slower but sometimes
  1895.             faster than modes 2-3.
  1896.         4,5,8,9 ==> left-to-right digit generation.
  1897.         6-9 ==> don't try fast floating-point estimate
  1898.             (if applicable).
  1899.  
  1900.         Values of mode other than 0-9 are treated as mode 0.
  1901.  
  1902.         Sufficient space is allocated to the return value
  1903.         to hold the suppressed trailing zeros.
  1904.     */
  1905.  
  1906.     int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
  1907.         j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
  1908.         spec_case, try_quick;
  1909.     Long L;
  1910. #ifndef Sudden_Underflow
  1911.     int denorm;
  1912.     unsigned Long x;
  1913. #endif
  1914.     Bigint *b, *b1, *delta, *mlo, *mhi, *S;
  1915.     double d2, ds, eps;
  1916.     char *s, *s0;
  1917.     static Bigint *result;
  1918.     static int result_k;
  1919.  
  1920.     if (result) {
  1921.         result->k = result_k;
  1922.         result->maxwds = 1 << result_k;
  1923.         Bfree(result);
  1924.         result = 0;
  1925.         }
  1926.  
  1927.     if (word0(d) & Sign_bit) {
  1928.         /* set sign for everything, including 0's and NaNs */
  1929.         *sign = 1;
  1930.         word0(d) &= ~Sign_bit;    /* clear sign bit */
  1931.         }
  1932.     else
  1933.         *sign = 0;
  1934.  
  1935. #if defined(IEEE_Arith) + defined(VAX)
  1936. #ifdef IEEE_Arith
  1937.     if ((word0(d) & Exp_mask) == Exp_mask)
  1938. #else
  1939.     if (word0(d)  == 0x8000)
  1940. #endif
  1941.         {
  1942.         /* Infinity or NaN */
  1943.         *decpt = 9999;
  1944.         s =
  1945. #ifdef IEEE_Arith
  1946.             !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
  1947. #endif
  1948.                 "NaN";
  1949.         if (rve)
  1950.             *rve =
  1951. #ifdef IEEE_Arith
  1952.                 s[3] ? s + 8 :
  1953. #endif
  1954.                         s + 3;
  1955.         return s;
  1956.         }
  1957. #endif
  1958. #ifdef IBM
  1959.     d += 0; /* normalize */
  1960. #endif
  1961.     if (!d) {
  1962.         *decpt = 1;
  1963.         s = "0";
  1964.         if (rve)
  1965.             *rve = s + 1;
  1966.         return s;
  1967.         }
  1968.  
  1969.     b = d2b(d, &be, &bbits);
  1970. #ifdef Sudden_Underflow
  1971.     i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
  1972. #else
  1973.     if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) {
  1974. #endif
  1975.         d2 = d;
  1976.         word0(d2) &= Frac_mask1;
  1977.         word0(d2) |= Exp_11;
  1978. #ifdef IBM
  1979.         if (j = 11 - hi0bits(word0(d2) & Frac_mask))
  1980.             d2 /= 1 << j;
  1981. #endif
  1982.  
  1983.         /* log(x)    ~=~ log(1.5) + (x-1.5)/1.5
  1984.          * log10(x)     =  log(x) / log(10)
  1985.          *        ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
  1986.          * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
  1987.          *
  1988.          * This suggests computing an approximation k to log10(d) by
  1989.          *
  1990.          * k = (i - Bias)*0.301029995663981
  1991.          *    + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
  1992.          *
  1993.          * We want k to be too large rather than too small.
  1994.          * The error in the first-order Taylor series approximation
  1995.          * is in our favor, so we just round up the constant enough
  1996.          * to compensate for any error in the multiplication of
  1997.          * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
  1998.          * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
  1999.          * adding 1e-13 to the constant term more than suffices.
  2000.          * Hence we adjust the constant term to 0.1760912590558.
  2001.          * (We could get a more accurate k by invoking log10,
  2002.          *  but this is probably not worthwhile.)
  2003.          */
  2004.  
  2005.         i -= Bias;
  2006. #ifdef IBM
  2007.         i <<= 2;
  2008.         i += j;
  2009. #endif
  2010. #ifndef Sudden_Underflow
  2011.         denorm = 0;
  2012.         }
  2013.     else {
  2014.         /* d is denormalized */
  2015.  
  2016.         i = bbits + be + (Bias + (P-1) - 1);
  2017.         x = i > 32  ? word0(d) << 64 - i | word1(d) >> i - 32
  2018.                 : word1(d) << 32 - i;
  2019.         d2 = x;
  2020.         word0(d2) -= 31*Exp_msk1; /* adjust exponent */
  2021.         i -= (Bias + (P-1) - 1) + 1;
  2022.         denorm = 1;
  2023.         }
  2024. #endif
  2025.     ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
  2026.     k = (int)ds;
  2027.     if (ds < 0. && ds != k)
  2028.         k--;    /* want k = floor(ds) */
  2029.     k_check = 1;
  2030.     if (k >= 0 && k <= Ten_pmax) {
  2031.         if (d < tens[k])
  2032.             k--;
  2033.         k_check = 0;
  2034.         }
  2035.     j = bbits - i - 1;
  2036.     if (j >= 0) {
  2037.         b2 = 0;
  2038.         s2 = j;
  2039.         }
  2040.     else {
  2041.         b2 = -j;
  2042.         s2 = 0;
  2043.         }
  2044.     if (k >= 0) {
  2045.         b5 = 0;
  2046.         s5 = k;
  2047.         s2 += k;
  2048.         }
  2049.     else {
  2050.         b2 -= k;
  2051.         b5 = -k;
  2052.         s5 = 0;
  2053.         }
  2054.     if (mode < 0 || mode > 9)
  2055.         mode = 0;
  2056.     try_quick = 1;
  2057.     if (mode > 5) {
  2058.         mode -= 4;
  2059.         try_quick = 0;
  2060.         }
  2061.     leftright = 1;
  2062.     switch(mode) {
  2063.         case 0:
  2064.         case 1:
  2065.             ilim = ilim1 = -1;
  2066.             i = 18;
  2067.             ndigits = 0;
  2068.             break;
  2069.         case 2:
  2070.             leftright = 0;
  2071.             /* no break */
  2072.         case 4:
  2073.             if (ndigits <= 0)
  2074.                 ndigits = 1;
  2075.             ilim = ilim1 = i = ndigits;
  2076.             break;
  2077.         case 3:
  2078.             leftright = 0;
  2079.             /* no break */
  2080.         case 5:
  2081.             i = ndigits + k + 1;
  2082.             ilim = i;
  2083.             ilim1 = i - 1;
  2084.             if (i <= 0)
  2085.                 i = 1;
  2086.         }
  2087.     j = sizeof(unsigned Long);
  2088.     for(result_k = 0; sizeof(Bigint) - sizeof(unsigned Long) + j <= i;
  2089.         j <<= 1) result_k++;
  2090.     result = Balloc(result_k);
  2091.     s = s0 = (char *)result;
  2092.  
  2093.     if (ilim >= 0 && ilim <= Quick_max && try_quick) {
  2094.  
  2095.         /* Try to get by with floating-point arithmetic. */
  2096.  
  2097.         i = 0;
  2098.         d2 = d;
  2099.         k0 = k;
  2100.         ilim0 = ilim;
  2101.         ieps = 2; /* conservative */
  2102.         if (k > 0) {
  2103.             ds = tens[k&0xf];
  2104.             j = k >> 4;
  2105.             if (j & Bletch) {
  2106.                 /* prevent overflows */
  2107.                 j &= Bletch - 1;
  2108.                 d /= bigtens[n_bigtens-1];
  2109.                 ieps++;
  2110.                 }
  2111.             for(; j; j >>= 1, i++)
  2112.                 if (j & 1) {
  2113.                     ieps++;
  2114.                     ds *= bigtens[i];
  2115.                     }
  2116.             d /= ds;
  2117.             }
  2118.         else if (j1 = -k) {
  2119.             d *= tens[j1 & 0xf];
  2120.             for(j = j1 >> 4; j; j >>= 1, i++)
  2121.                 if (j & 1) {
  2122.                     ieps++;
  2123.                     d *= bigtens[i];
  2124.                     }
  2125.             }
  2126.         if (k_check && d < 1. && ilim > 0) {
  2127.             if (ilim1 <= 0)
  2128.                 goto fast_failed;
  2129.             ilim = ilim1;
  2130.             k--;
  2131.             d *= 10.;
  2132.             ieps++;
  2133.             }
  2134.         eps = ieps*d + 7.;
  2135.         word0(eps) -= (P-1)*Exp_msk1;
  2136.         if (ilim == 0) {
  2137.             S = mhi = 0;
  2138.             d -= 5.;
  2139.             if (d > eps)
  2140.                 goto one_digit;
  2141.             if (d < -eps)
  2142.                 goto no_digits;
  2143.             goto fast_failed;
  2144.             }
  2145. #ifndef No_leftright
  2146.         if (leftright) {
  2147.             /* Use Steele & White method of only
  2148.              * generating digits needed.
  2149.              */
  2150.             eps = 0.5/tens[ilim-1] - eps;
  2151.             for(i = 0;;) {
  2152.                 L = d;
  2153.                 d -= L;
  2154.                 *s++ = '0' + (int)L;
  2155.                 if (d < eps)
  2156.                     goto ret1;
  2157.                 if (1. - d < eps)
  2158.                     goto bump_up;
  2159.                 if (++i >= ilim)
  2160.                     break;
  2161.                 eps *= 10.;
  2162.                 d *= 10.;
  2163.                 }
  2164.             }
  2165.         else {
  2166. #endif
  2167.             /* Generate ilim digits, then fix them up. */
  2168.             eps *= tens[ilim-1];
  2169.             for(i = 1;; i++, d *= 10.) {
  2170.                 L = d;
  2171.                 d -= L;
  2172.                 *s++ = '0' + (int)L;
  2173.                 if (i == ilim) {
  2174.                     if (d > 0.5 + eps)
  2175.                         goto bump_up;
  2176.                     else if (d < 0.5 - eps) {
  2177.                         while(*--s == '0');
  2178.                         s++;
  2179.                         goto ret1;
  2180.                         }
  2181.                     break;
  2182.                     }
  2183.                 }
  2184. #ifndef No_leftright
  2185.             }
  2186. #endif
  2187.  fast_failed:
  2188.         s = s0;
  2189.         d = d2;
  2190.         k = k0;
  2191.         ilim = ilim0;
  2192.         }
  2193.  
  2194.     /* Do we have a "small" integer? */
  2195.  
  2196.     if (be >= 0 && k <= Int_max) {
  2197.         /* Yes. */
  2198.         ds = tens[k];
  2199.         if (ndigits < 0 && ilim <= 0) {
  2200.             S = mhi = 0;
  2201.             if (ilim < 0 || d <= 5*ds)
  2202.                 goto no_digits;
  2203.             goto one_digit;
  2204.             }
  2205.         for(i = 1;; i++) {
  2206.             L = d / ds;
  2207.             d -= L*ds;
  2208. #ifdef Check_FLT_ROUNDS
  2209.             /* If FLT_ROUNDS == 2, L will usually be high by 1 */
  2210.             if (d < 0) {
  2211.                 L--;
  2212.                 d += ds;
  2213.                 }
  2214. #endif
  2215.             *s++ = '0' + (int)L;
  2216.             if (i == ilim) {
  2217.                 d += d;
  2218.                 if (d > ds || d == ds && L & 1) {
  2219.  bump_up:
  2220.                     while(*--s == '9')
  2221.                         if (s == s0) {
  2222.                             k++;
  2223.                             *s = '0';
  2224.                             break;
  2225.                             }
  2226.                     ++*s++;
  2227.                     }
  2228.                 break;
  2229.                 }
  2230.             if (!(d *= 10.))
  2231.                 break;
  2232.             }
  2233.         goto ret1;
  2234.         }
  2235.  
  2236.     m2 = b2;
  2237.     m5 = b5;
  2238.     mhi = mlo = 0;
  2239.     if (leftright) {
  2240.         if (mode < 2) {
  2241.             i =
  2242. #ifndef Sudden_Underflow
  2243.                 denorm ? be + (Bias + (P-1) - 1 + 1) :
  2244. #endif
  2245. #ifdef IBM
  2246.                 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
  2247. #else
  2248.                 1 + P - bbits;
  2249. #endif
  2250.             }
  2251.         else {
  2252.             j = ilim - 1;
  2253.             if (m5 >= j)
  2254.                 m5 -= j;
  2255.             else {
  2256.                 s5 += j -= m5;
  2257.                 b5 += j;
  2258.                 m5 = 0;
  2259.                 }
  2260.             if ((i = ilim) < 0) {
  2261.                 m2 -= i;
  2262.                 i = 0;
  2263.                 }
  2264.             }
  2265.         b2 += i;
  2266.         s2 += i;
  2267.         mhi = i2b(1);
  2268.         }
  2269.     if (m2 > 0 && s2 > 0) {
  2270.         i = m2 < s2 ? m2 : s2;
  2271.         b2 -= i;
  2272.         m2 -= i;
  2273.         s2 -= i;
  2274.         }
  2275.     if (b5 > 0) {
  2276.         if (leftright) {
  2277.             if (m5 > 0) {
  2278.                 mhi = pow5mult(mhi, m5);
  2279.                 b1 = mult(mhi, b);
  2280.                 Bfree(b);
  2281.                 b = b1;
  2282.                 }
  2283.             if (j = b5 - m5)
  2284.                 b = pow5mult(b, j);
  2285.             }
  2286.         else
  2287.             b = pow5mult(b, b5);
  2288.         }
  2289.     S = i2b(1);
  2290.     if (s5 > 0)
  2291.         S = pow5mult(S, s5);
  2292.  
  2293.     /* Check for special case that d is a normalized power of 2. */
  2294.  
  2295.     if (mode < 2) {
  2296.         if (!word1(d) && !(word0(d) & Bndry_mask)
  2297. #ifndef Sudden_Underflow
  2298.          && word0(d) & Exp_mask
  2299. #endif
  2300.                 ) {
  2301.             /* The special case */
  2302.             b2 += Log2P;
  2303.             s2 += Log2P;
  2304.             spec_case = 1;
  2305.             }
  2306.         else
  2307.             spec_case = 0;
  2308.         }
  2309.  
  2310.     /* Arrange for convenient computation of quotients:
  2311.      * shift left if necessary so divisor has 4 leading 0 bits.
  2312.      *
  2313.      * Perhaps we should just compute leading 28 bits of S once
  2314.      * and for all and pass them and a shift to quorem, so it
  2315.      * can do shifts and ors to compute the numerator for q.
  2316.      */
  2317. #ifdef Pack_32
  2318.     if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
  2319.         i = 32 - i;
  2320. #else
  2321.     if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
  2322.         i = 16 - i;
  2323. #endif
  2324.     if (i > 4) {
  2325.         i -= 4;
  2326.         b2 += i;
  2327.         m2 += i;
  2328.         s2 += i;
  2329.         }
  2330.     else if (i < 4) {
  2331.         i += 28;
  2332.         b2 += i;
  2333.         m2 += i;
  2334.         s2 += i;
  2335.         }
  2336.     if (b2 > 0)
  2337.         b = lshift(b, b2);
  2338.     if (s2 > 0)
  2339.         S = lshift(S, s2);
  2340.     if (k_check) {
  2341.         if (cmp(b,S) < 0) {
  2342.             k--;
  2343.             b = multadd(b, 10, 0);    /* we botched the k estimate */
  2344.             if (leftright)
  2345.                 mhi = multadd(mhi, 10, 0);
  2346.             ilim = ilim1;
  2347.             }
  2348.         }
  2349.     if (ilim <= 0 && mode > 2) {
  2350.         if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
  2351.             /* no digits, fcvt style */
  2352.  no_digits:
  2353.             k = -1 - ndigits;
  2354.             goto ret;
  2355.             }
  2356.  one_digit:
  2357.         *s++ = '1';
  2358.         k++;
  2359.         goto ret;
  2360.         }
  2361.     if (leftright) {
  2362.         if (m2 > 0)
  2363.             mhi = lshift(mhi, m2);
  2364.  
  2365.         /* Compute mlo -- check for special case
  2366.          * that d is a normalized power of 2.
  2367.          */
  2368.  
  2369.         mlo = mhi;
  2370.         if (spec_case) {
  2371.             mhi = Balloc(mhi->k);
  2372.             Bcopy(mhi, mlo);
  2373.             mhi = lshift(mhi, Log2P);
  2374.             }
  2375.  
  2376.         for(i = 1;;i++) {
  2377.             dig = quorem(b,S) + '0';
  2378.             /* Do we yet have the shortest decimal string
  2379.              * that will round to d?
  2380.              */
  2381.             j = cmp(b, mlo);
  2382.             delta = diff(S, mhi);
  2383.             j1 = delta->sign ? 1 : cmp(b, delta);
  2384.             Bfree(delta);
  2385. #ifndef ROUND_BIASED
  2386.             if (j1 == 0 && !mode && !(word1(d) & 1)) {
  2387.                 if (dig == '9')
  2388.                     goto round_9_up;
  2389.                 if (j > 0)
  2390.                     dig++;
  2391.                 *s++ = dig;
  2392.                 goto ret;
  2393.                 }
  2394. #endif
  2395.             if (j < 0 || j == 0 && !mode
  2396. #ifndef ROUND_BIASED
  2397.                             && !(word1(d) & 1)
  2398. #endif
  2399.                     ) {
  2400.                 if (j1 > 0) {
  2401.                     b = lshift(b, 1);
  2402.                     j1 = cmp(b, S);
  2403.                     if ((j1 > 0 || j1 == 0 && dig & 1)
  2404.                     && dig++ == '9')
  2405.                         goto round_9_up;
  2406.                     }
  2407.                 *s++ = dig;
  2408.                 goto ret;
  2409.                 }
  2410.             if (j1 > 0) {
  2411.                 if (dig == '9') { /* possible if i == 1 */
  2412.  round_9_up:
  2413.                     *s++ = '9';
  2414.                     goto roundoff;
  2415.                     }
  2416.                 *s++ = dig + 1;
  2417.                 goto ret;
  2418.                 }
  2419.             *s++ = dig;
  2420.             if (i == ilim)
  2421.                 break;
  2422.             b = multadd(b, 10, 0);
  2423.             if (mlo == mhi)
  2424.                 mlo = mhi = multadd(mhi, 10, 0);
  2425.             else {
  2426.                 mlo = multadd(mlo, 10, 0);
  2427.                 mhi = multadd(mhi, 10, 0);
  2428.                 }
  2429.             }
  2430.         }
  2431.     else
  2432.         for(i = 1;; i++) {
  2433.             *s++ = dig = quorem(b,S) + '0';
  2434.             if (i >= ilim)
  2435.                 break;
  2436.             b = multadd(b, 10, 0);
  2437.             }
  2438.  
  2439.     /* Round off last digit */
  2440.  
  2441.     b = lshift(b, 1);
  2442.     j = cmp(b, S);
  2443.     if (j > 0 || j == 0 && dig & 1) {
  2444.  roundoff:
  2445.         while(*--s == '9')
  2446.             if (s == s0) {
  2447.                 k++;
  2448.                 *s++ = '1';
  2449.                 goto ret;
  2450.                 }
  2451.         ++*s++;
  2452.         }
  2453.     else {
  2454.         while(*--s == '0');
  2455.         s++;
  2456.         }
  2457.  ret:
  2458.     Bfree(S);
  2459.     if (mhi) {
  2460.         if (mlo && mlo != mhi)
  2461.             Bfree(mlo);
  2462.         Bfree(mhi);
  2463.         }
  2464.  ret1:
  2465.     Bfree(b);
  2466.     *s = 0;
  2467.     *decpt = k + 1;
  2468.     if (rve)
  2469.         *rve = s;
  2470.     return s0;
  2471.     }
  2472. #ifdef __cplusplus
  2473. }
  2474. #endif
  2475.